// eee.Sheffield.PZ.Imaging
//
// Copyright ?Ping Zou, 2007
// sg71.cherub@gmail.com

using System;
using System.Collections.Generic;
using System.Text;
using System.IO;
using eee.Sheffield.PZ.Math;

namespace eee.Sheffield.PZ.Imaging.ImageFilter
{
    // Zhou et al 1994
    public class ZhouMatchedFilter
    {
        public PZMath_vector MatchedFilter(PZMath_vector sl, double A, double squareSigma)
        {
            int width = sl.Size;
            double r = 1.96 * System.Math.Sqrt(squareSigma);
            int int_r = (int)r < (width - 1) / 2 ? (int)r : (width - 1) / 2 - 1;
            // output width
            int outputWidth = width - int_r > 0 ? width - int_r : 0;
            PZMath_vector output = new PZMath_vector(outputWidth);

            int lowerBound = int_r < outputWidth ? int_r : 0;
            int upperBound = outputWidth - 1;

            //if (lowerBound >= upperBound)
            //    throw new ApplicationException("lowerBound should be less than upperBound");

            for (int i = lowerBound; i <= upperBound; i++)
            {
                double iterm1 = 0.0;
                for (int j = 0; j < width; j++)
                {
                    double f_l = A * System.Math.Exp(-1.0 * (i - j) * (i - j) / 2.0 / squareSigma);
                    iterm1 += sl[j] * f_l;
                }
                output[i] = iterm1;
            }

            return output;
        } // ZhouMatchedFilter()       

        // find 2 gaussian function tails at +/- sigma, and r = 1.96 * sigma
        public double EstimateRadius(PZMath_vector sl, int c)
        {
            // total intensity
            double totalIntensity = 0.0;
            int width = sl.Size;
            for (int i = 0; i < width; i++)
                totalIntensity += sl[i];

            // start from centre, find +sigma position (right)
            int mean = c;
            int rightPos = mean;
            double tempSumIntensity = 0.0;
            double thresholdSumIntensity = totalIntensity * 0.341;  // 34.1%
            bool findRightPos = true;
            while (tempSumIntensity < thresholdSumIntensity)
            {
                tempSumIntensity += sl[rightPos];
                rightPos++;
                if (rightPos >= width)
                {
                    findRightPos = false;
                    break;
                }
            }

            // find -sigma position (right)
            int leftPos = mean;
            tempSumIntensity = 0.0;
            bool findLeftPos = true;
            while (tempSumIntensity < thresholdSumIntensity)
            {
                tempSumIntensity += sl[leftPos];
                leftPos--;
                if (leftPos < 0)
                {
                    findLeftPos = false;
                    break;
                }
            }
            // average
            double sigma = 0.0;
            double count = 0.0;
            if (findLeftPos)
            {
                sigma += (double)(mean - leftPos);
                count++;
            }
            if (findRightPos)
            {
                sigma += (double)(rightPos - mean);
                count++;
            }
            sigma /= count;

            // r = 1.96 * sigma
            double r = 1.96 * sigma;
            return r;
        } // EstimateRadius()

        public void Zhou(PZMath_vector sl,
            out double radius, out double centre, out double leftEdge, out double rightEdge)
        {
            int zhou_ln = sl.MaxIndex();         // ln
            int width = sl.Size;
            if (zhou_ln == 0)
                zhou_ln = 1;
            if (zhou_ln == width - 1)
                zhou_ln = width - 2;
            double zhou_A = sl[zhou_ln];          // A = g(ln)
            int zhou_lnM1 = zhou_ln - 1;                            // ln-1
            int zhou_lnA1 = zhou_ln + 1;                            // ln+1
            double zhou_R = System.Math.Log(sl[zhou_ln] / sl[zhou_lnA1]) /
                System.Math.Log(sl[zhou_lnM1] / sl[zhou_ln]);
            // R
            double zhou_C = (zhou_lnA1 * zhou_lnA1 - zhou_ln * zhou_ln -
                zhou_R * (zhou_ln * zhou_ln - zhou_lnM1 * zhou_lnM1)) /
                2 / (zhou_R * (zhou_lnM1 - zhou_ln) - (zhou_ln - zhou_lnA1));
            // C
            double zhou_squareSigma = ((zhou_ln - zhou_C) * (zhou_ln - zhou_C) -
                (zhou_lnM1 - zhou_C) * (zhou_lnM1 - zhou_C)) / 2 /
                System.Math.Log(sl[zhou_lnM1] / sl[zhou_ln]);
            // sigma ^ 2
            if (zhou_squareSigma < 0.0)
                zhou_squareSigma = 0.0;

            PZMath_vector zhou_output = MatchedFilter(sl, zhou_A, zhou_squareSigma);
            int zhou_m = zhou_output.MaxIndex();
            double zhou_std = sl.StdDev;
            //double zhou_r = 1.96 * zhou_std;
            double zhou_r = EstimateRadius(sl, zhou_m);

            PZMath_vector zhou_edgePoints = new PZMath_vector(2);
            zhou_edgePoints[0] = zhou_m - zhou_r;
            zhou_edgePoints[1] = zhou_m + zhou_r;

            // output
            radius = zhou_r;
            centre = (double)zhou_m;
            leftEdge = zhou_edgePoints[0];
            rightEdge = zhou_edgePoints[1];
        } // Zhou()        

        public void ZhouRecord(PZMath_vector sl, string fileName)
        {
            FileStream stream = new FileStream(fileName, FileMode.Create, FileAccess.Write);
            StreamWriter writer = new StreamWriter(stream);


            int zhou_ln = sl.MaxIndex();   // ln
            int width = sl.Size;
            if (zhou_ln == 0)
                zhou_ln = 1;
            if (zhou_ln == width - 1)
                zhou_ln = width - 2;
            double zhou_A = sl[zhou_ln];          // A = g(ln)
            writer.WriteLine(zhou_A);

            int zhou_lnM1 = zhou_ln - 1;                            // ln-1
            int zhou_lnA1 = zhou_ln + 1;                            // ln+1
            double zhou_R = System.Math.Log(sl[zhou_ln] / sl[zhou_lnA1]) /
                System.Math.Log(sl[zhou_lnM1] / sl[zhou_ln]);
            // R
            double zhou_C = (zhou_lnA1 * zhou_lnA1 - zhou_ln * zhou_ln -
                zhou_R * (zhou_ln * zhou_ln - zhou_lnM1 * zhou_lnM1)) /
                2 / (zhou_R * (zhou_lnM1 - zhou_ln) - (zhou_ln - zhou_lnA1));
            writer.WriteLine(zhou_C);

            // C
            double zhou_squareSigma = ((zhou_ln - zhou_C) * (zhou_ln - zhou_C) -
                (zhou_lnM1 - zhou_C) * (zhou_lnM1 - zhou_C)) / 2 /
                System.Math.Log(sl[zhou_lnM1] / sl[zhou_ln]);
            // sigma ^ 2
            if (zhou_squareSigma < 0.0)
                zhou_squareSigma = 0.0;

            writer.WriteLine(zhou_squareSigma);


            // draw zhou's gaussian
            PZMath_vector zhouGaussian = new PZMath_vector(width);
            for (int i = 0; i < width; i++)
            {
                zhouGaussian[i] = zhou_A * System.Math.Exp(-1.0 *
                    (i - zhou_C) * (i - zhou_C) / 2 / zhou_squareSigma);
                writer.Write(String.Format("{0, -10:f}", zhouGaussian[i]));
            }
            writer.WriteLine();

            PZMath_vector zhou_output = MatchedFilter(sl, zhou_A, zhou_squareSigma);
            // writer zhou output
            int outputWidth = zhou_output.Size;
            for (int i = 0; i < outputWidth; i++)
            {
                writer.Write(String.Format("{0, -10:f}", zhou_output[i]));
            }
            writer.WriteLine();

            int zhou_m = zhou_output.MaxIndex();
            writer.WriteLine(zhou_m);
            double zhou_std = sl.StdDev;
            //double zhou_r = 1.96 * zhou_std;
            double zhou_r = 1.96 * System.Math.Sqrt(zhou_squareSigma);
            writer.WriteLine(zhou_r);
            PZMath_vector zhou_edgePoints = new PZMath_vector(2);
            zhou_edgePoints[0] = zhou_m - zhou_r;
            zhou_edgePoints[1] = zhou_m + zhou_r;

            writer.Flush();
            writer.Close();
            stream.Close();
        } // ZhouRecord()     
    }
}
